Double Life of Methanol: Experimental Studies and Nonequilibrium Molecular-Dynamics Simulation of Methanol Effects on Methane-Hydrate Nucleation

We have investigated systematically and statistically methanol-concentration effects on methane-hydrate nucleation using both experiment and restrained molecular-dynamics simulation, employing simple observables to achieve an initially homogeneous methane-supersaturated solution particularly favorable for nucleation realization in reasonable simulation times. We observe the pronounced “bifurcated” character of the nucleation rate upon methanol concentration in both experiments and simulation, with promotion at low concentrations and switching to industrially familiar inhibition at higher concentrations. Higher methanol concentrations suppress hydrate growth by in-lattice methanol incorporation, resulting in the formation of “defects”, increasing the energy of the nucleus. At low concentrations, on the contrary, the detrimental effect of defects is more than compensated for by the beneficial contribution of CH3 in easing methane incorporation in the cages or replacing it altogether.

statistical mechanics ensemble associated to some value of a suitable set of observables, which can be the same or different from the ones one follows along the time evolution of the system.
The system can be prepared in an initial macroscopic condition by, for example, restrained molecular dynamics, ReMD, [3][4][5][6] which is described in the next section.
More in detail, let us denote a generic observable which is a function of the phase space (Γ) Γ. For the time-dependent case, the observable is an ensemble average in phase space of ( ) ( Γ) (1) bracket and L(t) the Liouville operator. Since numerical evaluation of f (Γ, t) cannot be computed explicitly, apart from simple cases, we apply the Onsager-Kubo relation where S (t) is the time-evolution operator, S*(t) is its adjoint, and f (Γ, t 0 ) is the PDF at the initial time t 0 . In so doing, we can compute the time-dependent average as the ensemble ( ) average over the initial PDF f (Γ, t 0 ) of the observable O(Γ(t)) at the point Γ (t), corresponding to the evolution in time of the initial phase-space point Γ (t 0 ). In practice, in the D-NEMD approach, we sample initial positions and velocities from the initial PDF f (Γ, t 0 ) by ReMD, where k is a tunable parameter. The corresponding canonical PDF is Noting that , so that in the limit of , we observe the right-hand side of eqn. S3 going to eqn. S1, and we obtain k   . Restraining molecular dynamics must not be confused with molecular dynamics using to relax the initial condition; In ReMD, using the terminology of statistical mechanics, time is used for sampling, while in the relaxation trajectory time preserves its genuine meaning.
We wish to sample nucleation events at different supersaturated conditions of methane and (relatively low) solvated-methanol concentrations; thus, we need to introduce observables counting number of solvated methane and methanol molecules in water solution, n CH4 and n CH3OH , respectively. The estimator n CH4 (Γ) is defined as where is a characteristic function associated to the a th sub cell. In other words, the molecule a  contributes with 1 to n CH4 when it is in the a th sub cell and 0 otherwise. This definition follows that given by Irving and Kirkwood for a microscopic field. 7 In our case, we define a unique sub cell, , as a central orthogonal sub volume delimited along the z-axis by two perpendicular flat  boundaries, located at z 1 and z 2 with z 1 < z 2 , while length and width are the same of the simulation cell. For this cell, a simple form of the characteristic function is where z i is the component of the i th particle along the z axis in the ordinary space, and is ( )   the Heaviside step function. It is worth stressing that the Heaviside step function is a discontinuous function of its domain, resulting to impulsive forces on atoms crossing the boundary of the cell due to the biasing term. This problem is usually addressed by replacing  by a smooth analytical approximation, here where λ is a parameter controlling the smoothness of the approximation. 8,9 S5 At high methane-supersaturation conditions, we note that the solution tends not only to relax towards the final hydrate state, but a significant part of solvated methane is involved in the nucleation of pure methane bubbles. This is mainly due to the low solubility of methane in water. The consequences are essentially twofold: (1) the homogeneity of the solution is not preserved, and (2) the concentration of solvated methane decreases in solution. Hence, we need to introduce a further parameter to avoid the formation of pure methane bubbles in the sub cell . Following previous successful works, [10][11][12] we exploit the difference in the coordination  number of the pair CH 4 -CH 4 at distance equal to 5.5 Å (in the following text denoted by CN CH4−CH4 ) between the bulk liquid solution, and the supercritical fluid of pure methane in the bubble. 10,13 With neighboring cutoff distance 5.5 Å, CN CH4−CH4 is close to zero for the bulk solution whilst it is greater in the pure methane bubble. Thus, CN CH4−CH4 can be used to discern/prevent the presence of pure methane bubbles in solution. Consistently with n CH4 , we define a suitable estimator of CN CH4−CH4 as where the argument of the summation is the smooth analytical approximation of already ( )   introduced in eqn. S8. R cut is equal to 5.5 Å, and N c is a suitable scalar factor. Using the ReMD approach, we restrain the value of CN CH4−CH4 around zero so that the homogeneity of the solution and the concentration of solvated methane are preserved during the time trajectory. It is worth stressing that the extended Hamiltonian introduced in eqn. S2 within the ReMD approach is now composed of three biasing quadratic potential terms, which restrain the time evolution around the target values of n CH4 , n CH3OH and CN CH4−CH4 , respectively.
To apply D-NEMD, we need to sample a proper stationary condition ensemble to assess the initial PDF f (Γ, t 0 ). In our case, the initial conditional ensemble describes a state consisting of an aqueous supersaturated phase of solvated methane. Thus, we need a suitable observable able to characterise the initial state and distinguish it from the final hydrate state of equilibrium.
Cluster-analysis hydrate-recognition methods provide non-analytical functions of phase space, and, consequently, cannot be used in the ReMD approach. Therefore, we use a revised analytical version of the parameter F 4 , already proposed by Rodger et al. 14 where N c is a scalar factor and is the dihedral angle, defined as the angle between two , , , i j k l  planes: the first plane is determined by the two non-collinear vectors and and the second  difference between the most (60°, 180°, and 300°) and least (120° and 240°) probable angles. 15 Thus, F 4 in bulk water is zero. In all crystalline structures where water molecules form almost planar 5-and 6-membered rings, F 4 adopts non-zero values, since the most probable angle is 180°. 14 To sample the initial PDF f (Γ, t 0 ) via D-NEMD, we add another biasing quadratic potential term to the extended Hamiltonian of eqn. S2 to restrain the parameter F 4 around zero, corresponding to the bulk solution. We stress that this term is active only during sampling of the initial PDF f (Γ, t 0 ).

Computational Details
We prepared 16 different computational samples arising from the combination of four concentrations of methane (χ CH4 =0.038, 0.044, 0.052, 0.058, denominate samples A, B, C and D) and four concentrations of methanol (χ CH3OH =0, 0.008, 0.016, 0.024) in water (see Table   S1). Samples are arranged in orthogonal simulation boxes, elongated along the z-axis, and are composed of two phases. The first phase is an aqueous solution of solvated methane and methanol at the given respective concentrations, and it is in the central part of the simulation box with respect to the z-axis. The second is a methane liquid phase located on the sides of the periodic simulation box with respect to the z-axis. The interfaces between the two phases are flat and perpendicular relative to the z-axis. For all simulation boxes, the water-methane interface was positioned at z 1 = −42 Å and z 2 = 42 Å. Further, we denominate the central orthogonal sub-volume between z φ1 = −32 Å and z φ2 = 32 Å on which we impose a prescribed concentration of methane and methanol. All samples were filled up randomly with the methane and water molecules, with numbers specified in Table S1 respecting   In this article we introduced a simple model for methanol, which is compatible and consistent with the water/methane model proposed by Jacobson and Molinero. 18 Our methanol molecule consists of the rigid union between a mW water unit, a point particle, and a united-atom CH 4 unit. The mW unit is a Stillinger-Weber point particle with a typical tetrahedral coordination with other mW units, modeling the typical tetrahedral coordination among water molecules in liquid water, ices and clathrates. The united-atom CH 4 unit interacts with mW water molecules S8 through the two-body repulsive term of the Stillinger-Weber potential; the parameters of this interaction have been tuned to reproduce some of the characteristics of methane/water mixtures. 18 Thus, bonding a CH 4 unit to an mW one prevents the formation of any hydrogen bonding with other mW water units along the mW-CH 4 direction, which reduces the maximum number of hydrogen bonds that can be formed by the methanol to three, like for real CH 3 OH.
In addition, this simplistic methanol model preserves the characteristic to be made of a hydrophilic head and a (short) hydrophobic tail. Summarizing, our mW-CH 4 methanol model, though very simple, presents the main characteristics of the real molecule.
To validate this model, we performed MD simulations of i) water/methanol solutions and ii) sI methane clathrate hydrate incorporating methanol in the structure. These calculations have been performed with both the models used in present work and with all-atoms force fields. For allatom calculations we used TIP3P and TIP4P water, 19 and OPLS 20,21 methane and methanol.
To validate the methanol model used in the present work we compare methanol oxygen-water oxygen (partial) pair correlation function g(r) in the liquid water/methanol mixture and sI methane clathrate incorporating methanol (Fig. S1). This allows us to validate the capability of our simple but effective CH 4 -mW methanol model to mimic the hydrogen bonding of all-atoms force field both in terms of structure and energetics. Indeed, it can be shown that the pair correlation function is related to the so-called potential of the mean force, 24 also known as the Landau free energy of the distance between two atoms/particles of the given species, i.e., the effective pair interaction between methanol oxygen and water oxygen in the given system:

Hydrate formation and mass-balance-based determination of hydrate-conversion yield
The experimental apparatus for hydrate formation employed a pressure vessel fabricated using 316 stainless steel with internal volume of approximately 340 cm 3 (Figs. S3 and S4). The vessel was agitated using a tilting shaker. A pressure transducer with an uncertainty of 0.2 MPa, was used to measure pressure, whilst a thermocouple with an accuracy of ±0.1 K was inserted into the cell to measure the inner temperature, with temperature/pressure readings every 2 s. Prior to each run, the vessel was cleaned thoroughly and sterilised by washing with ethanol solution (25 wt%). Each hydrate-formation experiment began with charging the equilibrium cell with approximately 20 cm 3 of deionised water. The main system was cooled to the desired temperature of 2.9 °C (to mimic seafloor temperatures), via the temperature-control system.

S10
Once at the desired temperature, the cell was evacuated for 3 min by vacuum pump to remove any residual air, and then pressurised to the desired pressure of 120 bar using pure methane.
Due to inevitable Joule-Thomson thermal contraction, the cell pressure was decreased slightly by decreasing the temperature; however, within less than 10 minutes, the temperature stabilises and remains so until the end of the chose hydrate-formation period. Then, a continuous slow pressure decline was observed during the hydrate crystal-growth stage (always under the constant-volume conditions), after nucleation. In practice, however, there is some small temperature fluctuation during hydrate formation (i.e., the thermal trace of hydrate formation) due to its exothermic nature, but this is countered continually by the temperature-control system. Here, the average temperature of the production régime (i.e., the temperature plateau) is considered as the starting temperature of the reaction. The system was kept at (or very near) the desired temperature of 3.5 °C for 24 hours to gauge the yield towards hydrate over this period.
The 'yield' for conversion to hydrate during formation is calculated based on the number of absorbed moles of gas into the liquid/solid phase, i.e., by monitoring gas-phase pressure drop continuously on a mass-balance basis. Naturally, the first step in this number-of-gas-phasemoles-from-pressure determination lies in defining accurate the de-facto compressibility factor of the methane in the system, measured on our own system and tailored specifically by its slight temperature gradients and thermal inertia, in terms of the readily measurable temperature and pressure data, as detailed in the Supporting Information. Of course, we also do make use of 'reference' literature values as a function of temperature and pressure. The number of gas moles in the chamber can then be extracted using available data via PV = zn RT.
Knowing the number of absorbed/released moles of methane, in addition to typical methanehydrate cage-occupancy levels of 90% in the present P/T range, 25 allows for the percentage yield to methane hydrate to be calculated.
Estimation of compressibility for the experimental determination of the conversion yield As mentioned above, 'yield' of conversion is calculated based on the number of absorbed moles of gas into the liquid/solid phase, which require determination of de-facto compressibility factor of the methane in the system. For this purpose, in five separate runs, the chamber was loaded with methane gas (and nothing else, after thorough cleaning of the type described above) at ~80 to 130 bar (spanning the entire gamut, and more, of recorded gas-phase pressure during hydrate-formation cycles) for up to five different pressure values. While all the inlet/outlet valve were closed, the reactor was heated up slowly from 3 to 5 °C whilst pressure and temperature S11 were recording by computer, and these data for the five points were selected to calculate the compressibility-factor value, z. Due to no completely and utterly definitive measurements of the number of loaded moles into the system, the working, effective z-value of one P/T pair was taken from ref. 26 as the correct one: the exact number of gas moles in the chamber can be extracted using this hypothesis, via PV = zn RT. Given that this number of moles, n, is the same for all the other P/T points, the de-facto z value of those other four points can be therefore inferred via the same equation over the five independent runs, and the mean and variance taken.
These sets of lowest-variance points with three coordinates (temperature, pressure and z-value) have been plotted as a 3D plot while the equation of a fitted surface to these points was obtained using a rational Taylor model, 27 as -(cf. Fig. S5). This surface-fit equation was employed later to calculate the effective, operational compressibility factor of methane at various pressures and temperatures in the hydrate-formation runs from mass balances on thus-inferred gas-phase-number-of-moles data (from the gas-phase pressure) after the first 10 minutes of thermal equilibration, taking into account the temperature-variation of methane absorption in liquid with literature data for Henry's-Law constants for methane. 28 S12